


********************************************************************************
********************************************************************************
** compare prior tests 
********************************************************************************
********************************************************************************

** 1. for each date t, get indicator for had test in [t-15, t-9]
	use health/test_daily_panel, clear
	merge m:1 sid using health/demographics, assert(2 3) keep(1 3) nogen ///
		keepusing(sid age_wide)
	sort sid date
	drop if age_wide == 0 | missing(age_wide)
		
	forvalues f = 9/15 {
		gen future`f' = date + `f'
	}
	list sid date in 1
	reshape long future, i(sid date) 
	
	format future %td
	list future tested if sid == 221
	keep sid future tested age_wide
	duplicates drop
	rename future date
	rename tested prior_test	
	save health/prior_tests, replace 

** 2. for each date t, get indicated for had test in [t-2, t+4]
	use health/test_daily_panel, clear
	merge m:1 sid using health/demographics, assert(2 3) keep(1 3) nogen ///
		keepusing(sid age_wide)
	sort sid date
	drop if age_wide == 0 | missing(age_wide)
		
	forvalues c = 1/7 {
		local f = `c'-3
		gen future`c' = date + `f'
	}
	list sid date in 1
	reshape long future, i(sid date) 
	
	format future %td
	list future tested if sid == 221
	keep sid future tested age_wide
	duplicates drop
	rename future date
	rename tested admit_tests
	save health/admit_tests, replace 


** 3. count up number of people tested in [t-2, t+4] and [t-15, t-9]
	
	use health/prior_tests
	merge 1:1 sid date using health/admit_tests, nogen 
	egen week = cut(date), at(`=mdy(3,13,2020)'(7)`=mdy(12,22,2020)')	
	drop if missing(week)
	collapse (max) prior_test admit_tests, by(sid week age_wide)	
	collapse (sum) prior_test, by(week age_wide)
	
	merge m:1 age_wide using data/age_wide_counts, assert(3) nogen
	
	gen prior_rate = prior_test / n
	collapse (mean) prior_rate [fw=n], by(week)
	
	
	gen prior_rate_pop = string(prior_rate, "%4.3f") + " & " 
	keep week prior_rate_pop
	save health/prior_rate_pop, replace 

** 4. merge the hospital data with prior test data	
	foreach group in any_icli no_icli any_major {
		use health/inp_clean if admit_date >= mdy(3,1,2020), clear
		
		if "`group'" == "no_icli" gen no_icli = any_icli == 0
		
		sort sid admit_date	
		by sid: keep if _n==1 & `group'==1
		
		merge m:1 sid using health/demographics, keep(1 3) assert(2 3) nogen ///
			keepusing(sid age_wide)
		drop if missing(age_wide)  | age_wide == 0
		
		rename admit_date date
		merge 1:1 sid date using health/prior_tests, keep(1 3) 
		replace prior_test = 0 if _merge == 1
		drop _merge
		
		merge 1:1 sid date using health/admit_tests, keep(1 3) 
		replace admit_test = 0 if _merge == 1
		
		sum test_tight admit_test if _merge ~= 2
		drop _merge 

		gen month = month(date)

		egen week = cut(date), at(`=mdy(3,13,2020)'(7)`=mdy(12,22,2020)')	
		drop if missing(week)
		
		collapse (mean) prior_test (count) n_test=prior_test, by(week age_wide)
		
		merge m:1 age_wide using data/age_wide_counts
		bysort week: egen tt = total(n)
		gen wt = n/tt
		gen avar = wt*prior_test*(1-prior_test)/n_test
		list age_wide prior_test n_test n wt if week == mdy(3,27,2020)				
		collapse (mean) prior_test (rawsum) n_test avar [fw=n], by(week )
		list prior_test if week == mdy(3,27,2020)
		
		
		gen se = sqrt(avar)
		
		gen upper = prior_test  + 1.96*se
		gen lower= prior_test - 1.96*se
		replace lower = 0 if lower < 0 
		
		gen prior_rate_`group' = string(prior_test, "%4.3f") + " " 
		gen ci_`group' = "(" + string(lower, "%4.3f") + ", " + ///
			string(upper, "%4.3f") + ") &"
		
		rename se se_`group' 
		rename prior_test prior_rate_n_`group'
		keep week prior_rate_`group' ci_`group' se_`group' prior_rate_n_`group'

		save health/prior_rate_`group', replace
	}

	
** 3. stitch together and outsheet
	use health/prior_rate_pop, clear
	foreach group in  no_icli any_major any_icli {
		merge 1:1 week using health/prior_rate_`group', assert(3) nogen
	}

	gen year = year(week)
	gen month = month(week)
	gen day = day(week)
	tostring year month day, replace
	replace day = "0" + day if length(day)==1
	gen weeks = month + "/" + day + "/" + year + " & " 
	replace ci_any_icli = subinstr(ci_any_icli, "&", "\\", .)
	
	list weeks prior_rate_pop prior_rate_no_icli ci_no_icli ///
		prior_rate_any_major ci_any_major ///
		prior_rate_any_icli ci_any_icli, noobs clean
	
	outsheet weeks prior_rate_pop prior_rate_no_icli ci_no_icli ///
		prior_rate_any_major ci_any_major ///
		prior_rate_any_icli ci_any_icli using tables/prior_rates.tex, ///
		replace noquote nonames delimit(" ")



** graphical version
replace prior_rate_pop = subinstr(prior_rate_pop, "&", "",.)
destring prior_rate* ,replace

gen upper_no_icli = prior_rate_no_icli + 1.96*se_no_icli
gen lower_no_icli = prior_rate_no_icli - 1.96*se_no_icli
replace lower_no_icli = 0 if lower_no_icli<0
gen upper_any_major = prior_rate_any_major + 1.96*se_no_icli
gen lower_any_major = prior_rate_any_major - 1.96*se_no_icli
replace lower_any_major = 0 if lower_any_major < 0

foreach t in no_icli any_major {
		
	if "`t'" == "no_icli" local tt "Non-ICLI hospitalizations"
	if "`t'" == "any_major" local tt "Clear cause hospitalizations"
		
	# delimit ;
	twoway
		(line prior_rate_pop week )		
		(line prior_rate_`t' week, color(maroon) )
		(rarea upper_`t' lower_`t' week, color(maroon%20))
		(pci .045 `=mdy(11,13,2020)' .061 `=mdy(11,13,2020)', color(navy))
		(pci .017 `=mdy(5,13,2020)' .042 `=mdy(5,13,2020)', color(maroon))
			,
		graphregion(color(white) margin(r+5))
		legend(off)
		xtitle("") ytitle("")
		xlabel(`=mdy(3,13,2020)' "3/13" `=mdy(4,13,2020)' "4/13"
			`=mdy(5,13,2020)' "5/13" `=mdy(6,13,2020)' "6/13"
			`=mdy(7,13,2020)' "7/13" `=mdy(8,13,2020)' "8/13" 
			`=mdy(9,13,2020)' "9/13" `=mdy(10,13,2020)' "10/13"
			`=mdy(11,13,2020)' "11/13" `=mdy(12,13,2020)' "12/13")
		name(`t', replace)
		title("`tt'")
		ylabel(0(0.02)0.08)
		text(.062 `=mdy(11,14,2020)' "Population", place(w) color(navy))
		text(.043 `=mdy(5,13,2020)' "Hospitalizations", place(n) color(maroon))
		
	;
	# delimit cr
}

graph combine no_icli any_major, cols(2) graphregion(color(white)) ///
	xsize(8) ysize(3) iscale(1)

graph export figures/prior_rates.pdf, replace 
	

